#include "DXUT.h"
#include "PCGGonioFilesIO.h"
#include "../Utility/D3D/D3DUtilityFunction.h"

CPCGGonioFilesIO::CPCGGonioFilesIO(void):m_iCurWavelengthBand(400)
{
}


CPCGGonioFilesIO::~CPCGGonioFilesIO(void)
{
}

void CPCGGonioFilesIO::ReadGonioFile(char* strFileName, unsigned int iBand)
{
	m_iCurWavelengthBand = iBand;
	
	FILE* fp = fopen(strFileName,"r");

	m_arrBRDFData.RemoveAll();
	while(!feof(fp))
	{
		PCGGonioMeasuredData data;
		float fThetaOut, fTheta2, fThetaIn, fBRDF;
		fscanf(fp,"%f %f %f %f\n",&fThetaOut,&fTheta2,&fThetaIn,&fBRDF);
		data.m_fThetaIn = fThetaIn; 
		data.m_fThetaOut = fThetaOut;
		data.m_fBRDF = fBRDF;

		m_arrBRDFData.Add(data);
	}
	fclose(fp);
}


bool CPCGGonioFilesIO::ComputeNDFInverse(ID3D11Device* pd3dDevice, TCHAR* strImgName, char* str1DPlotFile, unsigned int iWidth, unsigned int iHeight, float fRotateAngle)
{
	if(m_arrBRDFData.GetSize()==0)
		return false;

	float radianperdegree = D3DX_PI/180.0f;
	//loop all samples
	NDFCountElem* pNDFBuf = new NDFCountElem[iWidth*iHeight];
	//-90 to 90, 1 degree interval
	NDFCountElem* pNDFPlot = new NDFCountElem[180];

	for (unsigned int i = 0; i < m_arrBRDFData.GetSize(); i++)
	{
		//convert the theta, phi to 3d vector
		float inTheta = m_arrBRDFData[i].m_fThetaIn*radianperdegree;
		float inPhi = (0+fRotateAngle)*radianperdegree;
		float outTheta = m_arrBRDFData[i].m_fThetaOut*radianperdegree;
		float outPhi = (0+fRotateAngle)*radianperdegree;
			
		//printf("%f %f %f %f \n",m_AngleMappingArray[iCurIdx].fThetaIn,m_AngleMappingArray[iCurIdx].fPhiIn + m_OrgGonioDataArray[i].fRotateAngle,m_AngleMappingArray[iCurIdx].fThetaOut,m_AngleMappingArray[iCurIdx].fPhiOut + m_OrgGonioDataArray[i].fRotateAngle);
		//smith terms --> ignore the incident/receive theta > 60 degrees
		if(fabs(inTheta) > D3DX_PI/3.0f || fabs(outTheta) > D3DX_PI/3.0f)
			continue;
		//printf("%f %f %f %f \n",m_AngleMappingArray[iCurIdx].fThetaIn,m_AngleMappingArray[iCurIdx].fPhiIn + m_OrgGonioDataArray[i].fRotateAngle,m_AngleMappingArray[iCurIdx].fThetaOut,m_AngleMappingArray[iCurIdx].fPhiOut + m_OrgGonioDataArray[i].fRotateAngle);
		/*
		x= rcos(phi)sin(theta)
		y= rsin(phi)sin(theta)
		z = rcos(theta)
		*/
		D3DXVECTOR3 vInVec;
		vInVec.x = cos(inPhi)*sin(inTheta);
		vInVec.y = sin(inPhi)*sin(inTheta);
		vInVec.z = cos(inTheta);

		D3DXVECTOR3 vOutVec;
		vOutVec.x = cos(outPhi)*sin(outTheta);
		vOutVec.y = sin(outPhi)*sin(outTheta);
		vOutVec.z = cos(outTheta);

		D3DXVECTOR3 vHalfVec;
		D3DXVECTOR3 vTmp = 0.5f*(vInVec+vOutVec);
		D3DXVec3Normalize(&vHalfVec,&vTmp);

		//compute the theta for half vector
		float fThetaM = acosf(vHalfVec.z);
		if(fRotateAngle == 0)
		{
			if(vHalfVec.x < 0) 
				fThetaM *= -1;
		}
		else if(fRotateAngle == 90)
		{
			if(vHalfVec.y < 0) 
				fThetaM *= -1;
		}
			 
		//printf("%f \n", fThetaM/D3DX_PI*180);
		//convert to degree
		fThetaM = fThetaM/D3DX_PI*180; //convert to degree
		//printf("%f,", fThetaM);
		int idxArr = floor(fThetaM+0.5f) + 90 - 1; //index
		//printf("%d \n", idxArr);
		pNDFPlot[idxArr].iCount++;

		//compute the fresenl term
		//http://www.filmetrics.com/refractive-index-database/Stainless+Steel
		//Refractive Index of Stainless Steel
		//For a typical sample of Stainless Steel the refractive index and extinction coefficient at 632.8 nm are 2.757 and 3.792.
		//float fFresnel = CD3DUtilityFunction::ComputeFresnelTermUsingSchlick(2.757, vInVec, vHalfVec);
		float fFresnel = D3DXFresnelTerm(cos(inTheta),2.757);
		//printf("fFresnel=%f,fFresnel=%f\n",fFresnel,fFresnel1);
		//Geometry term == 1.0 --> smith

		//assume the surface sample perpendicular to the zenith
		D3DXVECTOR3 vNorm(0,0,1);

		//                   Fresnel(l,h) * Geometry(l,v,h) * NDF(h)
		// Microfacet BRDF = ------------------------------------------
		//                    4*(dot(n,l))*(dot(n,v))
		// let's take the wavelength 640 nm
		float brdf = m_arrBRDFData[i].m_fBRDF;
		//float brdf = m_OrgGonioDataArray[i].RfRadiance390nm;
		//float brdf = m_OrgGonioDataArray[i].RfRadiance730nm;
		float ndotl = D3DXVec3Dot(&vNorm,&vInVec);
		float ndotv = D3DXVec3Dot(&vNorm,&vOutVec);
		float ndf = brdf*4*ndotl*ndotv/fFresnel;
		//printf("%f \n", ndf);
		pNDFPlot[idxArr].fVal += ndf;

		if (ndf < 0)
		{
			printf("ERROR: negative ndf!\n");
			return false;
		}

		//find the right position in the buffer for half vector
		vHalfVec.x = vHalfVec.x/2.0f + 0.5f; //(-1,1) --> (0,1); 
		vHalfVec.y = vHalfVec.y/2.0f + 0.5f;
		int iPosX = floorf(vHalfVec.x * iWidth);
		int iPosY = floorf(vHalfVec.y * iHeight);
		//printf("%d %d\n",iPosX,iPosY);
		unsigned int idx = iPosY*iWidth + iPosX;

		pNDFBuf[idx].iCount++;
		pNDFBuf[idx].fVal += ndf;
	}

	//create texture for saving
	float fMaxVal = DBL_MIN;
	for (int j = 0; j < iHeight; j++)
		for(int i = 0; i < iWidth; i++)
		{
			//printf("%f\t", pAccumBuf[j*iWidth+i]);
			//average NDF values
			if(pNDFBuf[j*iWidth+i].iCount>0)
				{
					pNDFBuf[j*iWidth+i].fVal /= pNDFBuf[j*iWidth+i].iCount;
					printf("%f,%d\n",pNDFBuf[j*iWidth+i].fVal, pNDFBuf[j*iWidth+i].iCount);
				}

			if(pNDFBuf[j*iWidth+i].fVal > fMaxVal)
				fMaxVal = pNDFBuf[j*iWidth+i].fVal;
		}

		HRESULT hr = S_OK;
		// Create texture
		D3D11_TEXTURE2D_DESC TexDesc;
		TexDesc.Width              = iWidth;
		TexDesc.Height             = iHeight;
		TexDesc.MipLevels          = 1;
		TexDesc.ArraySize          = 1;
		TexDesc.Format             = DXGI_FORMAT_R8G8B8A8_UNORM;

		DXGI_SAMPLE_DESC RealSampleDesc;
		RealSampleDesc.Count = 1;
		RealSampleDesc.Quality = 0;
		TexDesc.SampleDesc         = RealSampleDesc;

		TexDesc.Usage              = D3D11_USAGE_DYNAMIC;
		TexDesc.BindFlags          = D3D11_BIND_SHADER_RESOURCE;
		TexDesc.CPUAccessFlags     = D3D11_CPU_ACCESS_WRITE;
		// If they request mipmap levels, it's nice to be able to autogenerate them.
		TexDesc.MiscFlags          = 0;

		ID3D11Texture2D* pTex;
		V(pd3dDevice->CreateTexture2D(&TexDesc, NULL, &pTex));
		// Update the description with the read number of mipmaps, etc.
		pTex->GetDesc(&TexDesc);
		//
		// map for writing
		ID3D11DeviceContext* pd3dImmediateContext;
		pd3dDevice->GetImmediateContext(&pd3dImmediateContext);
		D3D11_MAPPED_SUBRESOURCE MappedResource; 
		V(pd3dImmediateContext->Map(pTex, 0, D3D11_MAP_WRITE_DISCARD, 0, &MappedResource));

		RGBA8Bytes* pPixel = (RGBA8Bytes*)(MappedResource.pData);
		for (int j = 0; j < iHeight; j++)
			for(int i = 0; i < iWidth; i++)
			{
				unsigned int idx = j*iWidth + i;
				if(pNDFBuf[idx].fVal < 0.000001f)//==0
				{
					(pPixel + idx)->r = 0;
					(pPixel + idx)->g = 0;
					(pPixel + idx)->b = 0;
					(pPixel + idx)->a = 255;
				}
				else{
					D3DXVECTOR3 color = CD3DUtilityFunction::GetColour(pNDFBuf[idx].fVal,0,fMaxVal);
					(pPixel + idx)->r = color.x * 255;
					(pPixel + idx)->g = color.y * 255;
					(pPixel + idx)->b = color.z * 255;
					(pPixel + idx)->a = 255;
				}

			}
			pd3dImmediateContext->Unmap(pTex,0);
			V(D3DX11SaveTextureToFile(pd3dImmediateContext,pTex,D3DX11_IFF_BMP,strImgName));

			// finish writing	
			SAFE_RELEASE(pTex);
			SAFE_RELEASE(pd3dImmediateContext);
			SAFE_DELETE_ARRAY(pNDFBuf);

			//save plot files
			float sumNDF = 0;
			for (int i = 0; i < 180; i++)
			{
				if(pNDFPlot[i].iCount > 0)
					pNDFPlot[i].fVal /= pNDFPlot[i].iCount;

				sumNDF += pNDFPlot[i].fVal;
			}
				
			FILE* fp = fopen(str1DPlotFile,"w");
			for (int i = 0; i < 180; i++)
			{
				if(pNDFPlot[i].fVal > 0)
				{
					int newidx = (i+1) - 90;
					fprintf(fp,"%d %f\n",newidx,(pNDFPlot[i].fVal/sumNDF)); //normalize ndf into [0,1]
				}
			}
			fclose(fp);
			SAFE_DELETE_ARRAY(pNDFPlot);
}